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ABSTRACT 

We present the systemic Console, a new all-in-one, general-purpose software package for the analysis 
and combined multiparameter fitting of Dopplcr radial velocity (RV) and transit timing observations. 
We give an overview of the computational algorithms implemented in the Console, and describe the 
tools offered for streamlining the characterization of planetary systems. We illustrate the capabilities 
of the package by analyzing an updated radial velocity data set for the HD128311 planetary system. 
HD128311 harbors a pair of planets that appear to be participating in a 2:1 mean motion resonance. 
We show that the dynamical configuration cannot be fully determined from the current data. We find 
that if a planetary system like HD128311 is found to undergo transits, then self-consistent Newtonian 
fits to combined radial velocity data and a small number of timing measurements of transit midpoints 
can provide an immediate and vastly improved characterization of the planet's dynamical state. 
Subject headings: Extrasolar Planets, Data Analysis and Techniques 



1. INTRODUCTION 

During the past decade, the characterization of ex- 
trasolar planets has become a major branch of Astron- 
omy. The field is energized by a variety of ground 
and space-based detection programs that are meeting 
with increasing success. In the past year, the cen- 
sus of extrasolar planets has exceeded 300, and plan- 
ets have now been successfully detected using a va- 
riety of_te£hmo 1 ues, in^^ velocity 
(e.g. iMavor fe Quelozlll995t lUdrv et al.ll2007bD. trans it 
photometry (e.g. iHenrv et al.l l2000t ICharbonneau et alJ 
120001 120071). microlensing (iBennettj I2009T) astrometrv 
(iBenedict etal . 2002; B ean fe Seifahrtll2009D . stellar pul- 



sations (ISilvotti et alJ 2007f) and even di rect imaging, 
(IChauvin et all 120051 ; iKalas et all 120081 ; iMarois et all 
12008ft . 

The radial velocity method has been used to discover 
more than 75% of the known planets, and continues to 
be a dominant techni que, both in terms of its contin- 
ued productivity (e.g. iFischer et all l2005f ) and its abil- 
ity to accurately probe planetary architect ures into the 
vicin i ty of the terrestrial mass region (e. 
2001 lLovis et all [200l lUdrv et all [2007 
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20091 ). A number of planets that were initially detected 
using radial velocity (e.g. HD 209458b, HD 189733b, HD 
149026b, Gl 436b, HD17156b and HD80606b) have been 
later shown to transit as a result of follow-up photometry, 
and because the parent stars of these planets arc bright, 
follow-up characterizations with a variety of method s 
have been extremely valuable fe.g. lDeming et aLll2005| ). 

The planets that have been detected with the radial ve- 
locity technique comprise a complicated and non-uniform 
sample. Some systems such as U psilon Andromedae 
dButler et all HflOflL 12001 . GJ 876 (jMarcv et all fl998l . 
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20061 ) have had multiple planets subject to very accurate 
orbital characterization within uniform, well-sampled 
data sets. Other sys tems, for example Epsilon Eridani 
([Benedict et al.ll2006h . draw their support from a variety 
of observational sources and in some cases have orbital 
parameters that are significantly uncertain. Indeed, it 
is difficult to draw a firm boundary between detections 
that are secure, and those that may be subject to serious 
revision or even elimination. 

In addition to the large amount of observational work 
that has gone into the detection of extrasolar planets, 
there is a parallel effort by theorists to explain the emerg- 
ing distributions of planets within the context of theories 
of planetary formation and evolution. This work spans 
a wide variety of bases, but a unifying principle is that 
much of it depends on the raw data being supplied by 
the catalog of extrasolar planets, and therein lies a dif- 
ficulty. Dynamicists have traditionally dealt with plane- 
tary orbital elements that are known to exquisite preci- 
sion. As far back as the Eighteenth century, the orbital 
elements of the solar system planets were known with an 
accuracy well in excess of our current orbital determina- 
tions for extrasolar planets. Theoretical interpretations 
of the extrasolar planetary data is sometimes made with- 
out full account of the highly varying signal-to-noise of 
the datasets that make up the catalog. This problem is 
exacerbated by the fact that there exists no continuously 
up-to-date compendium of known extrasolar planets in 
which all of the fits are derived using the same toolset 
of routines. The systemic collaboration has been estab- 
lished as an effort to solve this problem. 

The plan for this paper is as follows. In £j2j we describe 
the systemic Console. In SJ2 we show some sample appli- 
cations of the tools that are incorporated in the Console, 
with a particular emphasis on the planetary system or- 
biting HD128311 (|Vogt et all 120051 ). We show that our 
current radial velocity data set for this system is insuffi- 
cient for characterizing the resonant relation between the 
planets, and we demonstrate, using synthetic datasets, 
how the inclusion of transit timing data (were transits 
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to be detected) would almost immediately eliminate this 
degeneracy. As another example of the versatility of the 
code, we describe in Appendix [B] an automated pipeline 
(the systemic "backend" ) which runs on top of the same 
program to create a web application that analyses data 
sets and aggregates fits. In [J4j we describe the direc- 
tion of possible future work with the tools that we have 
developed, and conclude. 

2. THE SYSTEMIC CONSOLE 

The systemic Console is a downloadable software pack- 
age 4 that provides an intuitive graphical user interface 
for the fitting of planetary signatures, and an associated 
suite of dynamical analysis tools (Table [lj. It can also be 
used as a specialized, programmable calculator and run 
scripts in non-interactive mode to access its library of 
numerical routines. The program is written in the Java 
programming language for cross-platform portability. 

2.1. Radial Velocities 

The systemic Console allows for a choice between two 
modeling schemes. For the majority of the known extra- 
solar planetary systems, the planets do not experience 
significant dynamical interactions during the time range 
spanned by a set of radial velocity observations. In these 
cases, the radial velocity variation of the star can be rep- 
resented as a sum of N Keplerian orbits, each described 
by orbital elements (period P, mass A4, eccentricity e, 
mean anomaly M, and argument of pcriastron 137.) 

Summed Keplerians provide an adequate model for 
nearly all of the planetary systems that have been discov- 
ered to date (Appendix |A|) . Kepler's equation is rapidly 
solved using a simple iterative scheme, and hence mod- 
els can be quickly evaluated (see e.g. iFordl [200 9;. for a 
discussion of the current state-of-the-art). 

T here are, however, several exce ptions, notably GJ 
876 dRivera et al.|[2005h. HD202206 (ICorreia et al.l l2005l) 
and HD60532 ( jLaskar fc Correial l2009f) in which a self- 
consistent, or Newtonian fit is required. In these cases, 
planetary interactions are taken into account in the fit. 
and the Console adopts an N-planet model of the system 

^GA^(xi_Xj) ^ 

.7 = 1 



d 2 Xj 
dt 2 



with the integrations carried out using either 4th/5th or- 
der Runge-Kutta with adap tive timestep control or Her- 
mite 4th-order integration (|Press et al.l Tl992; Hu t et all 
Il995f h When an integrated model is adopted, a system is 
defined by the osculating orbital elements of the planets 
at the epoch of the first observation expressed in Jacobi 
coordinates (see Lee & Peale 2002). The user also has 
the option of providing an integration routine. 

Finally, the Console allows the velocity offsets between 
different data sources to be additional free parameters; 
this allows sources with different zero-point offsets (e.g. 
radial velocity surveys using different templates) to be 
combined in the fitting procedure. 

The Console carries out parameter minimization of the 
so-called reduced chi-square statistic 

N 
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of a fit; in the above expression, TV is the number of 
radial velocity data points, and M/a is the number of 
activated parameters, a\ . . . o,m- As a rule of thumb, a 
reduced Chi-square value near unity is indicative of a 
"good" fit to the data, suggesting that the model is a 
reasonable explanation of the data within the observa- 
tional errors. Typically, larger values usually signal an 
insufficient modeling of the data, whereas smaller values 
imply that the data has been over-fit. However, this rule 
is not exact, and should hence be applied with caution. 

2.2. Transits 

A rapidly growing number of planets (58 as of writing) 
with a favorably inclined orbital plane are being further 
characterized with transit timing data 5 . Transits enable 
direct estimations of planetary masses, radii and mean 
densities, tog ether with period and ph ase of the tran- 
siting planet (|Charbonneau et al.l I2007D . Considerable 
current interest is focused on detection of transit tim- 
ing variations (TTVs) which can point to the presence 
of additional perturbing bodies in a given system. 

When supplied to the Console, transits data (central 
primary and secondary transits timing) is included with 
the RV data in the following way. The Console searches 
for the best-fit orbital parameters by minimizing over the 
joint x 2 statistic 
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where Xrv represents the goodness-of-fit for the radial 
velocity component of the model, as described above, and 
Xt r is representative of the transit component. Ideally, 
one would fit together all of the radial velocity and tran- 
sit photometry data with a single model to jointly invert 
for the parameters that describe all available data. In 
the future, these capabilities will be incorporated into 
the Console. Much progress can still be made, however, 
by restricting our analysis to observed times of central 
transit with error-bars obtained from separate light curve 
analyses. These transit time data can then act as sep- 
arate constraints on the observed behavior of the sys- 
tem. To ease implementation, we compare the predicted 
and observed location of the planet at the observed time 
of central transit, rather than comparing transit times. 
Since the orbital velocities are not changing significantly 
with respect to the duration of the eclipse, the difference 
between these approaches is negligible. We thus use the 
following equation to define the goodness-of-fit statistic 
for the transit component of the model: 
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N 

E 
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where Sxi is the predicted separation perpendicular to 
the line of sight at the observed central transits ti, such 
that 

Sxi = \x*(ti) - x P (ti)\, i=l-N (5) 

The error on dxi is estimated from the error on ti as 
&8x.i = Vx,p&ti ■ While we do not explore it here, it is 
important to recognize that regularization of the fit may 
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5 Gary, B. 2009; http : //brucegary .net/AXA/x. html accessed 
13 March 2009 
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TABLE 1 
List of tools 



Name 



Menu command 



Description 



Fitting 

Levenberg-Marquardt 
Simulated Annealing 
De-trend 
Transit fitting 

Periodograms 

Lomb-Scarglc Pcriodogram 



Lomb-Scargle Pcriodogram of residuals View 



Edit -» Polish 

Edit — > Simulated Annealing 
Edit — » Detrend 
Options — » Transit fitting 



Periodogram of sampling 



Uncertainty estimation 

Bootstrap 



Markov Chain Monte Carlo 



F-test 

Dynamical analysis 

Dynamical evolution 

Transits prediction 



View — > Periodogram 

— > Periodogram of residuals 
View — » Periodogram of sampling 

View — » Bootstrap 

View — > Markov Chain Monte Carlo 

View — > F-test / F-test significance 



View 



View 



Orbital evolution and stability 
Transits prediction 



Multidimensional local optimization J2.3.2|. 
Multidimensional global optimization l|2.3,3[ l. 
Removes linear trends from the radial velocity data. 
Adds transits to the \ 2 statistics. 



Identifies periodicities in the full RV datasct 
and estimates FAPs. 

Identifies periodicities in the residual RV datasct 
and estimates FAPs. 

Identifies spurious periodicity peaks associated 
with uneven sampling of the radial velocities. 



Estimates uncertainties using the bootstrap routine; 
plot and ex port m arginal distributions of orbital 
parameters 112.4. 111 . 

Estimates uncertainties using the MCMC routine; 
check chain convergence; plot and e xport marginal 
distributions of orbital parameters H2.4.2I I. 



Tracks the fully integrated evolution of the orbital 
elements and the stability of the system. 
Calculates the distribution of central transit times 
for a given observational window. 



be warranted in this type of analysis (|Press "et~all[l99l) . 

6 

Since it is routinely possible to achieve small error 
bars on the central primary transits (100s for ground- 
based observations down to 10s for HST observations), a 
best fit found by the Console that includes transit timing 
may yield extremely precise determinations of the period 
and mean anomaly at epoch of the transit ing planet (e.g. 
IWittenmver et al.ll2005l ; iBean et al.|[200l 

De tection of central secondary eclipses (jDeming et al.l 
|2007| ) also places tight bounds on the eccentricity and ar- 
gument of pcriastron of the planet. This additional con- 
straint can break degeneracies present when RVs alone 
are used; for instance, it can discriminate between ec- 
centric single-planet systems and t wo-planet systems in a 
2:1 re sonance with circular orbits (|Anglada-Escude et al.l 
EH). 

Further afield, it can be possible to measure transit 
timing variations (TTV) in a dynamically interacting 
planetary configuration and infer t he orbital elements 
of a p erturbing, non-transiting body dHolman fc Murray! 

3n|l2 



of a p ( 
I2005H 



Aeol et all 2005; A gol fe Steffenll2007ft . 



2.3.1. 



2.3. Best-fit model estimation 
Periodograms and False Alarm Probabilities 



The Lomb-Scargle (LS) periodogram is an al gorithm 
for time series analysis of unevenly spaced data (jScargld 

6 Regularization is a formal statistical method of compromising 
between two distinct sources of information. This is accomplished 
by adding a relative weighting factor A in front of one of the com- 
ponents of the overall \ metric, where the value of A determines 
the relative importance of the two components of the goodness-of- 
fit. There are many different methods that can be used to choose 
an appropriate value for the weighting factor. In this work, we 
have implicitly chosen the value A = 1, corresponding to an equal 
weighting. 



fl98l iHorne fc Bahmmsl fl98l iPress et all [l992h . The 

LS periodogram is useful for rapidly identifying peri- 
odic signals in the observed data, and to residuals to 
a given fit, without having to fit for the other orbital 
parameters. The formula for an error-weighted peri- 
odogram P x (oj) as implemen ted in the Console is given in 
iGilliland fc Baliunasl ( 1987t ); the individual weights are 
taken to be Wj = 1/cJ- 

An advantage of this method is that its statistical prop- 
erties are well known and are conducive to the defini- 
tion of an analytic false alarm probability (FAP) associ- 
ated with each periodic signal. When the periodogram 
is normalized by the total variance Pq(lo) = P x {u)/o- , 
the estimated probability that a peak as high or higher 
would occur by chance is given by Pr(po, Nf) = 1 — [1 — 
cxp(— po)] Nf 7 where Nf is the effective number of fre- 
quencies. 

Finally, since the unequal spacing of the data can be 
a source of spurious periodicities (e.g. those associated 
with the synodic lunar month or yearly observational 
schedules), the Consol e also support s plotting of the 
power spectral window (Dccmina) ll975[ ) overlaid over the 
standard (non-error weighted) periodogram. 

2.3.2. Levenberg-Marquardt (local minimization) 

Given the observations and associated errors, the goal 
is to obtain a model configuration ybf (a 5N vector of 
orbital parameters) such that x 2 ( v b/) = mm yX 2 j this 
is usually reported as the "best-fit" solution. Typically, 
the Lomb-Scarglc pcriodogram is used to comb through 
periodicities in the data; periodicities are removed in or- 
der of decreasing half-amplitude K and optimized us- 
ing line-minimization. This procedure leads to a set of 
orbital parameters yo which is a rough approximation 
to the best-fit solution, and can be improved with si- 
multaneous multiparameter minimization. For a discus- 
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sion of the intricacies o f the Keplcrian fitting process, see 
ICumming et~aTl ()2008f ). 

Multidimensional parameter minimization can be car- 
ried out using the Levenberg-Marquardt algorithm (LM; 
iPress et al.lll992l ). Given the initial guess yo, the LM 
algorithm can quickly converge to a local minimum y'. 
Good convergence of the LM algorithm is conditional on 
the choice of the initial guess and a favorable geometry of 
the x 2 (y) surface: in particular, the algorithm is sensi- 
tive to rugged x 2 surfaces and can be prone to converging 
to non-optimal minima. 

2.3.3. Simulated Annealing (global minimization) 

So-called "global" minimization techniques attempt to 
avoid getting trapped in local minima by adding a de- 
gree of randomness at each iteration step, although at a 
muc h greater computa tional cost. Simulated annealing 
(SA; IPress et al.l 11992] ), by analogy to several thermo- 
dynamic processes in nature, defines an "energy" E as 
the objective function to minimize and allows for tem- 
perature fluctuations between states at different energies 
as dictated by the current temperature T n ; the tempera- 
ture T n is lowered with a (problem-dependent) scheduler. 
This algorithm is particularly appropriate for rugged x 2 
surfaces, or when the initial guess is sufficiently distant 
from the best-fit solution. 

In our problem, the objective function is clearly x 2 ( v )- 
Given a state y„, the algorithm selects a new configu- 
ration y rl +i; the new configuration is accepted and kept 
with a probability Pin — * n + 1) ~ cxp(— AE/T n ) if 
E n+ i > E n , and is always accepted if E n+ \ < E n (a 
downhill step). The temperature is subsequently up- 
dated according to the input scheduler, and the process is 
repeated until a target number of steps N is reached. The 
fact that uphill steps are sometimes accepted (according 
to the current temperature) lets the algorithm explore a 
larger portion of the parameter space and makes it less 
likely to get stuck in a narrow local minimum. The trial 
configuration y„+l is selected using a proposal distribu- 
tion, which is an easy-to-evaluate generator of trial con- 
figurations that picks a new set of parameters given the 
current set of parameters. The default function is a mul- 
tivariate Gaussian distribution centered on the current 
step y„; the variance fS^ can be chosen independently 
for each parameter. 

The algorithm requires that the following are config- 
ured from the user: 

1. temperature scheduler: the default scheduler de- 
creases T according to T n = T (l — n/N) a , where 
To and a are input parameters that dictate the ini- 
tial temperature and cooling rate. The optimal val- 
ues of To and a are problem-dependent and quite 
often may determine whether the routine success- 
fully recovers the true global minimum. 

2. generator of trial configurations: the default gen- 
erator is a Gaussian function centered around the 
current configuration, with the scale parameter 
vector (3^ given by the user (an initial value is sug- 
gested). 

Since the correct recovery of y&/ depends on appro- 
priate choices of Tg,a,N and (3^ that are not known 
a-priori, the Console allows several SA jobs to run in 



parallel, improving the chance of convergence to the best- 
fit model. Reconfigurations, in the form of occasionally 
jump-starting the routine with the best-ever solution, 
can also be beneficial to the success of the algorithm. 

Other global minimization schemes, s uch as 

Genetic Algorithms (e.g . ICharbonneaul 119951 : 
lLaughlin fc Chambers! f2001h . are being considered 
for inclusion in the Console's built-in array of tools. 
They can be easily implemented by the user using the 
routine library offered by the Console. 

Finally, we note that certain planetary systems such as 
HD80 606 (jLaughlin et all [20091 : lGilloril2009l : iPont et all 
|2009| ) include both photometric and spectroscopic data, 
and contain planets with high orbital eccentricities. In 
these cases, the connection between observable quanti- 
ties and the orbital and physical parameters is highly 
nonlinear, and a modeling framework that relies purely 
on x 2 minimization may have a difficult time recovering 
the correct system configuration. Future releases of the 
console will therefore incorporate the option of using a 
fully Baycsian approach to the fitting problem. 

2.4. Error estimation 

Radial- velocity searches are constantly pushing the en- 
velope towards lower and lower masses, frequently at the 
threshold of detection, with low signal-to-noise ratios. 
For this reason, once the best-fit parameters have been 
identified, it is vital to rigorously characterize their un- 
certainty. The Console offers two independent methods 
for estimating uncertainty: synthetic datasets refitting 
(bootstrap) and Markov Chain Monte Carlo (MCMC). 

2.4.1. Bootstrap 

The bootstrap procedure consists of drawing with re- 
placement from the observed data points (RV and central 
transits) and creating a number of synthetic data sets 
Af = i n- The Levenberg-Marquardt fitting procedure is 
then applied to each dataset, using the best-fit solution 
for the real dataset as the initial guess. The distribution 
of the obtained fitted parameters yf =1 ^ N yield an esti- 
mated a for the scatter of the orbital elements around 
the true intrinsic orbital parameters. 

Th e bootstrap algorithm is well known (|Press et al.l 
[1992^ and in common use for estimating planetary ele- 
ments uncertainties, but presents a number of disadvan- 
tages; namely, that it drives a local minimization routine 
(and is thus subject to the same pitfalls), and that it 
has a large computational burden. To partially improve 
on the first weakness, bootstrap can optionally be pre- 
ceded by a burn-in phase. The burn-in phase obtains a 
rough estimate of the scatter in the parameters by run- 
ning a short bootstrap phase. The error estimate is then 
used in the actual bootstrap run to perturb the best fit 
a set number of times (e.g. 10 times) per each synthetic 
dataset fitting; only the best-fitting final configuration 
is retained. This helps improving the reliability of the 
bootstrap routine in some cases. 

2.4.2. Markov Chain Monte Carlo 

Markov Cha in Monte Carlo (see, e.g., iFordl 120051 : 
lGregorvl2005bL for exoplanet related implementations) is 
an alternative method for estimating uncertainties that 
does not rely on minimization schemes. The MCMC 
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method generates a sequence (chain) of configurations 
yi that is sampled from the (unknown) probability dis- 
tribution /(y). The transition probability between two 
subsequent configurations y„ and y n +i is 



a(yn+i|y«) = min exp 



12 
Xn+1 



1 



(6) 



Assuming that the observational errors are accurately 
estimated and approximately Gaussian, this transition 
function assures that, after discarding an initial burn-in 
phase, the distribution of generated configurations will 
be sampled from the unknown probability distribution 
/• 

The algorithm consists of looping over the following 
steps, given an initial state yo: 

1. given a state y n and a Gaussian generator of trial 
states with scale parameters /3 M (see 12.3.5)) , gener- 
ate a trial state y'; 

2. accept the trial state y' with a probability a(y'|y„) 
and set y n +i = y', otherwise discard it (downhill 
steps are again always accepted); 

3. set ri = n + 1; 

until some convergence criterion of the chain is satis- 
fied. The MCMC algorithm guarantees convergence to 
the true distribution /(y), but can explore the parameter 
space inefficiently depending on the choice of or may 
not achieve satisfactory convergence within the chosen 
N steps. The convergence can be visually monitored by 
interactive plotting of the marginal distribution of the 
parameters. The acceptance rate of the MCMC algo- 
rithm can be interactively monitored as well; an optimal 
acceptance rate is ~ 0.25 (Gclman ct al. 2003). 

As with simulated annealing, multiple MCMC chains 
can be generated in parallel to provide comparison be- 
tween the results obtained with different choices of (3^ 
and chain length, which yield similar results within sta- 
tistical uncertainties if all chains have converged. More 
sophisticated B ayesian algorith ms, such as parallel tem- 
pering MCMC (lGregorvll2005ah . may be implemented by 
the user by exploiting the programmable interface of the 
Console. 

3. APPLICATIONS 

3.1. Resonance characterization in the HD128311 
system 

A high fraction of the detected extrasolar systems with 
multiple planet arc involved in near low-order mean mo- 
tion resonances (MMRs), with at least four of them 
(GJ876, HD82943, HD73526 and HD128311) being re- 
ported to engage in strong 2:1 resonances. Two planets 
are in a mean-motion resonance when the periods are in 
a ratio of small integers, and at least one of the reso- 
nant angles librates (i.e. it spans a range smaller than 
2tt). Resonant angles are linear combinations of w (ar- 
gument of periastron) and A = M + w (coplanarity is 
assumed). The relevant resonant angles for a 2:1 res- 
onance are 9\ = 2X 2 — X i — and Aw = w% — w\ 
(jMurrav fc Dermottl[2000T) . 

Radial veloci ty measurements for HD 1283 11 
(|Vogt et al.l |2005[ ) indicated that the system is locked 



in a 2:1 MMR, which ensures the long-term stability 
of the two giant planets. The best-fitting model was 
indefinitely stable, with the resonant argument B\ 
librating with an half-amplitude of about 60 degrees; 
a naive fit using Keplerian ellipses instead of the full 
N-body model is catastrophically unstable within about 
2,000 years. Orbital fits for the systems generated using 
a Monte Carlo procedure (similar to Section 12.4. ip 
yielded a proportion of about 60% stable systems with 
6*i librating and Azu circulating to about 40% with both 
arguments librating (apsidal co-rotation). The large 
stellar jitter (~ 9 m/s) and the relatively long periods 
of the two planets implies that models with different 
resonant configuration are equally likely given the radial 
velocity dataset. 

However, whether or not the system is in ap- 
sidal co-rotation is a crucial piece of information, 
since it can provide fundamental clues to the migra- 
tion and interaction history of the system. Scenar- 
ios of slow m igration and resonant capture into a 
2:1 M MR (e.g. iNelson fc Papaloizoul l2002t iLee fc Peald 
l2002t iBeauge et al.l |2006|) consistently result in sys- 
tems that are l ibrati ng in both resonant arguments. 
ISandor fc Kiev! (|2006f ). analyzing the specific case of 
HD128311, showed that after adiabatic migration and 
capture into MMR, the two planets are in apsidal co- 
rotation and have very small libration amplitudes. If 
a definitive prevalence of model fits not in apsidal co- 
rotation were ascertained, then the discrepancy has to 
be explained in terms of subsequent perturbative events 
(such as sudden termination of migration or planet- 
planet scattering) that happen after an adiabatic mi- 
gration pr ocess. Analogous studies have be en conducted 
for G J876 (jKlev et al.ll2005h and HD73526 (ISandor et all 
l2007h . 

It is therefore important to distinguish between the two 
resonant configurations (ideally, at the 90% confidence 
level or better); this requires a more precise determina- 
tion of the orbital parameters, which might be achieved, 
for instance, with additional RV measurements. 

For this purpose, we present a set of additional 
Doppler measurements taken between Jun e 2005 and 
May 2008 using the HIRES spectrometer (jVogt et all 
Il994h . Doppler measurements are taken using the stan- 
dard iodine cell technique (see iVoet et al.ll2009l for more 
details). Table [2] lists the updated Keck dataset, giving 
the time of each observation, the radial velocity and the 
internal uncertainties. 

3.2. Best fit 

We update the analysis of lVogt et al.1 (|2005h using the 
tools built in the Console for both the original data and 
the updated RVs presented in this paper. The Console 
is well suited to this task, since it can easily derive self- 
consistent fits (interactively) and do batch Monte-Carlo 
dynamical analyses on large sets of orbital parameters 
(non- interactively) . 

The two prominent periodicities in the IVogt et alJ 
(|2005h dataset were found using the integrated Lomb- 
Scargle periodogram. A self-consistent (Newtonian) 
best-fit was then derived using the Lcvenbcrg-Marquardt 
minimization routine; one of the built-in N-body inte- 
grators (Hermite) was used to derive the radial velocity 
curve for each choice of orbital parameters. The final 
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TABLE 2 

New Radial Velocities for HD128311 
(Sample: full table in electronic version) 



JD 


RV [m/s] 


Uncertainty [m/s] 


2450983.82690 


-12.95 


1.45 


2451200.13787 


-21.49 


1.92 


2451342.85836 


62.75 


2.05 


2451370.82904 


105.66 


1.88 


z4014Uy. i 4DDU 


1ZO. ( 1 


1 fiO 

1 .oz 


2451410.74909 


118.14 


2.01 


2451552.16457 


68.78 


1.85 


2451581.17009 


13.35 


1.64 


2451680.02544 


-60.10 


2.17 


2451974.16142 


62.03 


1.74 


2451982.15276 


32.30 


1.48 


2452003.02274 


12.76 


1.87 


2452003.90155 


29.10 


1.98 


2452005.13013 


27.90 


1.55 


2452061.87832 


-40.29 


1.54 


2452062.86745 


-11.26 


1.67 



TABLE 3 
Orbital fit parameters 





Fit A 


Fit B 




Fit C 


Period (d) 


466.6 [7.5] 


469.1 [3.3 




464.84 


909.5 [21.0] 


893.5 [6.2 




901.63 


Mass (Mj) 


1.59 [0.22] 


1.79 [0.17 




1.72 




3.19 [0.11] 


3.19 [0.08 




3.13 


Mean anomaly (deg) 


270.6 [31.9] 


282.2 [16.8 


] 


263.10 




192.0 [23.3] 


190.0 [13.7] 


193.33 


Eccentricity 


0.36 [0.07] 


0.33 [0.05] 


0.32 




0.20 [0.09] 


0.23 [0.05] 


0.20 


Long, of periastron (deg) 


73.8 [24.8] 


58.98 [19.6] 


78.04 




11.7 [20.0] 


4.54 [14.4 




6.59 



Note. — Fit A: integrated best-fit to the lVogt et al.l (12005ft Keck 
RV data. Fit B: integrated best-fit to the up dated RV data re- 
ported in this paper and the HET data reported infWittcnmyc r et al.| 
(2009). Fit C: orbital elements used to generate the synthetic 
datascts. All elements arc defined at epoch JD — 2450983.8269. 
Uncertainties arc reported in brackets. 

best-fit orbital parameters are listed as Fit A (Table \3§ . 
The uncertainties for each orbital parameter are found 
using the bootstrap routine on 10,000 synthetic dataset 
realizations. 

Subsequently, we derived the best-fit for the full up- 
dated Keck data (Table [2|) , together with the observa- 
tions taken with the Hobby-Eberly T elescope (HET) and 
reported in lWittenmver et al.l (|2009f ). The Lomb-Scargle 
periodogram and the associated analytic FAP estimates 
are shown in Figure [1] The Console can account for the 
zero-point offset and the velocity offset between the two 
datascts as two additional free parameters. The Newto- 
nian best-fit orbital parameters derived, however, result 
in a system that is unstable within 1000 years. There- 
fore, we generated a pool of alternative 5000 bootstrap- 
generated trial fits, checked each of them for stability 
within 10000 years and selected the best-fitting stable 
solution. Its orbital parameters and corresponding un- 
certainties are listed as Fit B (Table [3]). This model is 
protected by a 2:1 MMR, in which B\ librates with am- 
plitude ~ 60 deg and Azu circulates. The radial veloc- 
ity measurements and the star radial velocity curve are 
shown in Figure [2] 



180000 r- 
160000 - 
140000 - 
120000 - 

| 

g 100000 T 



80000 




1 10 100 1000 10000 

Days 



Fig. 1. — Lomb-Scargle periodogram for the combined Keck 
and HET datasets, as plotted by the Console. Analytical FAPs at 
levels 10 _1 (long dashed), 10 -2 (short dashed) and 10~ 3 (dotted) 
are overlaid. 
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-125.65 " 
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2450983.83 2451707.59 2452431.35 2453155.11 2453878.87 
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2450983.83 2451707.44 2452431.06 2453154.68 2453878.29 



Time 

Fig. 2. — Best-fit integrated solution to the RV data presented 
in this paper (blue) and the HET data (green) reported by Wit- 
tenmyer (orbital parameters listed as Fit B in Table \3^. 

3.3. Dynamical interactions 

Following the procedure detailed in lVogt et alj (|2005h . 
we took the two self-consistent two-planets fits (Fit A 
and Fit B) and applied a Monte-Carlo bootstrap proce- 
dure, in which new fits are derived by resampling with re- 
placements the radial velocity datasets. We created two 
Monte-Carlo generated libraries of 5000 self-consistent 
models for two radial velocities d atasets: the radial ve- 
locities listed in lVogt eT al. (2005]) and the updated Keck 
data reported in Table [2] For each of the two groups, 800 
fits, stable for at least 10 4 years 7 , were selected and inte- 
grated forward, recording the maximum eccentricity for 

7 For longer-term integrations, the builtin integration schemes 
(RK45 and 4th order Hermite) might not be sufficiently accurate 
and can be substituted by integration schemes supplied by the user. 
Alternati vely, the Console can be set up to drive packages such as 
SWIFT ( ht tp : //www . boulder . swri . edu/$\sim$hal/swif t . html ) 
or Mercury (Chambers & Migliorini 1997). 
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o i librating, Am circulating 
• | and Aru librating 
■ 6, and Am circulating 


2005 (Keck) 





: #JMf ' o 










■ 







2008 




Synthetic [102 RV] 
Synthetic [102 RV+ 4 transits] 
Synthetic [306 RV] 




e max ' irlner P lanet ) 

Fig. 3. — Maximum eccentricities observed during 10 4 yr inte- 
grations of self-consist e nt ob tained using the bootstrap routine for 
data from FVogt et all 1 120051) {top), data presented in this paper 
(middle) and synthetic data (bottom). Filled circles: scenarios in 
which both arguments librate. Open circles: scenarios in which 
9\ librates and Aro circulates. Gray squares: scenarios in which 
both arguments circulate. In the bottom panel, black and blue 
symbols are for models derived considering RV data only, whereas 
red symbols are for models considering RV and transits. 



TABLE 4 
Monte-Carlo analysis results 



Data 


R2 


Rl 


NR 


2005 (76 Keck RVs) 


281 [35%] 


489 [61%] 


30 


2008 (102 Keck RVs) 


160 [20%] 


618 [77%] 


22 


2008 (102 Keck + 78 HET RVs) 


180 22% 


615 [77%] 


5 


Fit C, 102 RVs 


603 [75%] 


197 [25%] 





Fit C, 102 RVs + 4 transits 


800 [100%] 








Fit C, 306 RVs 


743 [93%] 


52 [7%] 


5 



Note. — R2: resonant fits with both arguments librating. Rl: 
resonant fits with 8\ librating and Atu circulating. NR: fits have both 
arguments circulating. 

both planets and the amplitude of libration of both res- 
onant angles. The results of this analysis are shown in 
Figure S 

With the new radial velocity data, the percentage of 
model fits that arc stable and in apsidal co-rotation using 
the additional RVs falls slightly to 20%. A different run 
considering 1600 models also yields a similar percentage, 
confirming that the result is robust. The inclusion of the 
HET data also does not change our result significantly 
(Tabic [2]). Therefore, while we have strengthened the 
case for models of HD128311 that only librate in 9%, a 
secure determination of the libration amplitude of Azu 
might be obtained either by a transit monitoring cam- 
paign or yet additional RV measurements. 

3.4. Constraints by transits 
Although the a-priori geometric probability for transits 



Pi i 



P tr = 0.0045 



1AU\ fR* + R 



pi 



R, 



1 — e cos(-ij — w) 



1-e 2 



(7) 

(jBodenheimer et all l200l is very low for HD128311b 
(P tT fts 0.0032), given the high precision that can be 
achieved by the addition of transits to the % 2 budget, it 
is a worthwhile exercise as a proof of concept. Moreover, 
other resonant systems have higher transiting probabil- 
ities; for instance, planets GJ876b and c have a-priori 
transit probabilities *~ 1%, though the inclination of the 
system is unfavorable a nd no transits have been observed 
(iShankland et alJ[2006T l. 

We selected the best-fitting solution in apsidal corota- 
tion from the ensemble of systems generated by Monte- 
Carlo bootstrapping of the RVs presented in Table [2] (Fit 
C). The orbital elements are listed in Table [H Subse- 
quently, we created a synthetic dataset of RVs and tran- 
sits by integrating forward in time, using the N-body 
routines offered by the Console. The RVs are generated 
by sampling the radial velocity curve at the times listed 
in Table [5J the tabulated uncertainties and a jitter of 9 
m/s are added to the measurement. The central transit 
times dataset comprises four points, to which we added 
a Gaussian noise with amplitude 10~ 4 d (comparable to 
the uncertainties that ca n be achieved by g round-based 
transit observations; e.g. lAlonso et al.1120081 ). 

We repeated the analysis detailed in the previous sec- 
tion by bootstrapping exclusively the RV data (Table S]); 
this yields similar ratios, shifted to favor systems in ap- 
sidal corotation (similarly to the generating fit). 

As expected, the inclusion of the four central tran- 
sit times largely reduces the parameter space that can 
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be spanned by Monte-Carlo explorations. The large ex- 
cursions in x 2 an d the increased ruggedness of the x 2 
space makes the simple bootstrap algorithm, driving a 
Levenberg-Marquardt scheme, somewhat inefficient in 
fully exploring the allowed space of orbital parameters 
(as anticipated in Section f2 . 4 . 1 1) . We therefore used the 
Markov-Chain Monte Carlo routine supplied with the 
Console. A long chain of systems (5 x 10 5 ) was generated; 
the first 50000 systems were discarded and only 1 every 
fOO systems were retained, to minimize the correlation 
between subsequent chain elements. 

The tightness of the orbital parameter uncertainties 
thus generated (AP 1 /P 1 = 2.1 x 10~ 6 ; AP 2 /P 2 = 3 x 
10~ 6 d; AMi/Mi = 1.4 x 10~ 3 ; AM 2 /M 2 = 4.2 x 10~ 4 ; 
Aw i = 2.4 x 10~ 3 ; Azu 2 = 1.3 x 10~ 3 ) anticipates that 
the ratio of correctly recovered resonant configuration 
will be very high. In fact, with the addition of the four 
primary transits, all of the systems are correctly iden- 
tified in apsidal corotation (Table 2J. The maximum 
eccentricities achieved by the two planets (Figure [3|) are 
determined within 10 -3 . 

As a comparison, we ran the same procedure against 
204 additional RVs (a 30-year observation stretch), de- 
rived by sampling the integrated stellar radial velocity 
with the same schedule used for the Keck dataset. This 
large amount of additional RVs is required to identify the 
generating planetary system as apsidally corotating with 
a fraction >90% of models (Table HJ). 

4. DISCUSSION 

In this paper, we have described the features of the 
systemic software package. This software has been writ- 
ten with extensibility, portability and clarity as guiding 
principles, and is fully adequate for all but the most de- 
manding cxoplanet-related analysis tasks. The Console 
provides a uniform method for analyzing data stemming 
from a variety of sources (radial velocities surveys and 
transits) and allows the efficient recovery of the best- 
fitting stable planetary configuration, even in presence 
of strong mutual perturbations. It is provided for free to 
the scientific community. 

As an example application, we have analyzed an up- 
dated radial velocity dataset for the pair of resonat- 
ing pkiietshaibored by HD128311. As first noted 
by IVogt et all (|2005f ) , the orbital solution to this sys- 
tem is degenerate between apsidally corotating and non- 
apsidally corotating fits; the additional data sets do not 
break the degeneracy, owing to the large stellar jitter 
and long orbital periods. We have used an analysis of 
synthetic data sets to demonstrate that the detection of 
a transiting extrasolar planet system with planets par- 
ticipating in a low-order mean motion resonance, such 
as HD128311, would lead to a rapid determination of 



the libration widths of the resonant arguments and an 
attendant understanding in how such systems form and 
evolve. Additionally, our analysis shows that the pa- 
rameters of non-transiting planets can be very well con- 
strained through transit timing variations in presence of 
strong mutual interactions. As noted in Section (|2.2| . 
however, a more detailed analysis may be warranted (in 
particular regarding the issues of fit regularization and 
full photometry fitting) and will be the object of a follow- 
up paper. Finally, we showed that breaking the degen- 
eracy at a comparable level with radial velocities would 
require a prolonged observation campaign, of 30 years or 
more. 

We plan to improve the current feature set of the 
Console by (1) adding facilities for fully fitting the raw 
light curve data of a transit detection, (2) implementing 
more sophisticated routines for best-fit parameter and 
uncertainty estimation, and (3) allowing non-coplanar, 
inclined fits. We note that to date, nearly all of the plan- 
etary systems that have been detected with the Dopplcr 
radial velocity technique can be satisfactorily modeled 
(to the precision of the observations) using co-planar 
models with the inclinations assumed to be 90°. The 
Console's integration routines and internal system repre- 
sentations are fully three-dimensional, however, and so a 
forthcoming version will enable non-cop lanar fits and will 
accep t astrometric measurements (e.g. iBean fc Seifahrtl 
120091 ). With the advent of space missions such as SIM 
Lite and Gaia, there will be numerous opportunities to 
accurately discern the three-dimensional or bital configu- 
ratio ns of many nearby planetary systems (|Unwin et al.l 
l2008f) . Finally, signatures of less obvious effects in 
the spectroscopic and photometric d ata sets, such as 
those expected from general relativity (|Wu fc Goldreichl 
I2002T) or the excitati on of tidal modes in the host star 
(|Wu & Murra^l2003h . will require more sophisticated 
modelling to be properly taken into account. 



We would like to thank the participants in the sys- 
temic project for contributing a large amount of research 
effort toward the characterization of extrasolar planets. 
The results reported in this paper would not have been 
possible without their dedicated participation. We are 
grateful to Debra Fischer, Eric Ford, Man Hoi Lee, Doug 
Lin and Peter Jalowiczor for useful discussions, and the 
anonymous referee for a very thorough evaluation of the 
paper and several valuable suggestions. 

This research has been supported by the NSF through 
CAREER Grant AST-0449986, and by the NASA As- 
trobiology Institute through Grant NNG04GK19G. The 
Console software package may be downloaded for free at 
|http : / / www . oklo . org[ 



APPENDIX 

SUMMED KEPLERIANS MODEL 

When the perturbations between planets are negligible over the observational window, it is appropriate to model 
the radial velocity curve as a superposition of N Keplerian orbits of fixed orbital elements: 



N 

v r(t) = V"] Ki[cOs(vi + VJi) + Ci COS VJi] , 
i=l 



(Al) 
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where the radial velocity half-amplitude. Ki, of planet i is given by 

and where the true anomaly, Vi, is related to the eccentric anomaly, Ei, via 

(A3) 



tan 



1 + e, 
tan 

1 - e, 



Ei 
2 



The eccentric anomaly, Ei, in turn, can be expressed in terms of the mean anomaly Mi = 1~k / Pi(t — T per i,i) through 
Kepler's equation 

Mi = Ei- <* sin Ei (A4) 

(jMurrav fc Dermottl[2000f) . 

THE SYSTEMIC BACKEND 

The systemic backend is a web application that showcases the power of the Console as an automated engine for data 
analysis. It consists of a database of catalog information (stellar properties as well as RV and transit measurements) 
as published in the astronomical literature, and a catalog of model planetary fits for each star. For this purpose, it 
uses the Console as its main engine to perform a number of automatic data explorations, whereas the user-facing part 
uses standard "Web 2.0" tools (PHP, MySQL, Javascript and wikis) to present a coherent overview of the data. A 
public backend 8 is available as a proof-of-conccpt to foster collaboration within the broader community of exoplanet 
researchers and enthusiasts, and to present and maintain the catalog of fits to radial velocity and transit timing data 
for known planet-bearing stars. Each user has a personal data page and fit catalog, the possibility of commenting on 
other team member's fits, and can interact with other team members within a private and secure environment. A 
more powerful and customizable version is also available on request for use by planet hunter teams, and can be useful 
to maintain an integrated database of datasets and models in face of the increasing flux of RV and transit data. 

The fit catalog is scanned by a number of Console components, which continually sift through the uploaded fits 
in non-interactive mode. One component implements a bootstrap routine to calculate uncertainties on the orbital 
parameters of each fit; data from the bootstrap routine is stored in a database for creating scatter plots. Two other 
components check for dynamical instability over periods of 1,000 and 10,000 years, with stability defined by the rough 
criterion of requiring a smaller than 1% fractional change in semi- major axis with respect to the average semi- major axis 
observed during a full N-body integration. This step flags highly unstable planetary systems that experience ejections or 
collisions. Data from the integration is retained for plotting of orbital evolution and for future additional investigations. 
Dynamical data (orbital parameters, radial velocity data, fit parameters, stability, integrations, bootstrap results) is 
then transparently presented to the user as a set of web pages and can be aggregated and sliced using a web-based 
query system. 
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